library(broom)

#Model
model_plus1day <- "~ treatment_plus1day + age + ctr_gender + ctr_education +
  ctr_income + ctr_unemployed+ ctr_ethnicity + ctr_disability+
  daily_new_confirmed_per_million"

#DVs
dvs <- list(
  "likeCon" = "Like Cons (0-10)",
  "likeJohnson" = "Like Johnson (0-10)",
  "handleCorona" = "COVID-19 Performance (1-5)",
  "govtHandlelockdown" = "Lockdown Performance (1-5)",
  "convote" = "Vote Intention: Con (0-1)",
  "labvote" = "Vote Intention: Lab (0-1)"
)


parties <- c("1" = "Cons", "2" = "Lab", "10" = "Ind")

#Dataset to store the results
results_1day <- expand_grid(
  dv_name = names(dvs),
  party_id = names(parties)
) %>%
  mutate(estimate = NA_real_, std.error = NA_real_)

#Estimate the coefficients
for (i in seq_len(nrow(results_1day))) {
  
  dv  <- results_1day$dv_name[i]
  pid <- as.numeric(results_1day$party_id[i])
  
  df_tmp <- bes_w22_covid %>%
    filter(
      partyId == pid,
      date_diff >= -7,
      date_diff <= 6
    )
  
  if (dv %in% c("likeCon","likeJohnson","handleCorona","govtHandlelockdown")) {
    df_tmp <- df_tmp %>% filter(.data[[dv]] <= 10)
  }
  
  fml <- as.formula(paste(dv, model_plus1day))
  
  fit <- try(lm(fml, data = df_tmp), silent = TRUE)
  
  if (!inherits(fit, "try-error")) {
    est <- tidy(fit) %>% filter(term == "treatment_plus1day")
    if (nrow(est) == 1) {
      results_1day$estimate[i]  <- est$estimate
      results_1day$std.error[i] <- est$std.error
    }
  }
}

#Calculate CIs
df_scandal_plus1_robust <- results_1day %>%
  mutate(
    dv_name = dv_name,
    dv      = dvs[dv_name] |> unlist(),
    party   = parties[party_id],
    ul      = estimate + 1.96 * std.error,
    ll      = estimate - 1.96 * std.error,
    ul90    = estimate + 1.645 * std.error,
    ll90    = estimate - 1.645 * std.error,
    Party   = factor(party, levels = c("Cons","Lab","Ind")),
    DV      = factor(dv, levels = unlist(dvs))
  )

#Plot function
plot_plus1 <- function(df, dv_vars) {
  
  df %>%
    filter(dv_name %in% dv_vars) %>%
    ggplot() +
    geom_linerange(
      aes(x = Party, ymin = ll90, ymax = ul90),
      position = position_dodge2(width = 0.3),
      linewidth = 2,
      color = "black"
    ) +
    geom_pointrange(
      aes(x = Party, y = estimate, ymin = ll, ymax = ul),
      position = position_dodge2(width = 0.3),
      size = 2,
      linewidth = 1,
      color = "black"
    ) +
    facet_wrap(~DV) +
    geom_hline(yintercept = 0, linetype = 2) +
    labs(y = "Effect of Partygate Video") +
    theme_light() +
    theme(
      text = element_text(size = 30),
      strip.text = element_text(size = 26, color = "black"),
      strip.background = element_rect(fill = "grey90"),
      panel.grid = element_blank()
    )
}

#Figure J1
plus1_covid <- plot_plus1(
  df_scandal_plus1_robust,
  c("handleCorona", "govtHandlelockdown")
)
plus1_covid

#Figure J2
plus1_like <- plot_plus1(
  df_scandal_plus1_robust,
  c("likeCon", "likeJohnson")
)
plus1_like

#Figure J3
plus1_vote <- plot_plus1(
  df_scandal_plus1_robust,
  c("convote", "labvote")
)
plus1_vote

ggsave("plus1_like.eps",plus1_like,width = 14,height = 6,dpi = 1200)
ggsave("plus1_covid.eps",plus1_covid,width = 14,height = 6,dpi = 1200)
ggsave("plus1_vote.eps",plus1_vote,width = 14,height = 6,dpi = 1200)

